function [nIters, errRel] = testBlockMG(n, TOL, maxIt)
% A driver file to illustrate the usage of block-MG for Michael
%
% INPUT
%  n     : the number of cells, must be a power of 2.
%  TOL   : relative error threshold to stop the iteration.
%  maxIt : maximum number of iterations
%
% OUTPUT
%  nIters : the number of iterations performed
%  errRel : the relative error of the solution.
%
% example:
% [nIters, errRel] = testBlockMG(1024) should yield the following:
%  iter #1: max-norm of Residual=5.5435; max-norm of error=0.00085722
%  iter #2: max-norm of Residual=0.5113; max-norm of error=5.1334e-05
%  iter #3: max-norm of Residual=0.039407; max-norm of error=3.8525e-06
%  iter #4: max-norm of Residual=0.0029678; max-norm of error=2.8883e-07
%  iter #5: max-norm of Residual=0.00022272; max-norm of error=2.1654e-08
%  iter #6: max-norm of Residual=1.6702e-05; max-norm of error=1.6235e-09
%  iter #7: max-norm of Residual=1.2523e-06; max-norm of error=1.2173e-10
%  iter #8: max-norm of Residual=9.3891e-08; max-norm of error=9.1276e-12


% This subroutine assembles the Laplacian operator.
    function Lap = LaplacianOperator(n, Dim, h)
        % return the 2nd-order Laplacian operator
        %   with homogenuous Dirichlet BC.
        % assuming cell-centered variables and the vector is columnwise
        
        L   = sparse([],[],[],n,n);
        nm1 = ones(n-1,1);
        L   = L + diag(-2*ones(n,1), 0) + diag(nm1, 1) + diag(nm1, -1);
        % put in BC corrections.
        L(1,1)   = -3;
        L(n,n)   = -3;
        L = L/h/h;
        
        nP = n^Dim;
        if Dim==1
            Lap = L;
        elseif Dim==2
            In  = sparse(1:n, 1:n, ones(n,1), n, n);
            Lap = kron(L, In) + kron(In, L);
        else
            error('Unsupported value of Dim!');
        end
    end

if nargin<3
    maxIt = 30;
end
if nargin<2
    TOL = 1.e-6;
end

addpath('./core');
global mg_matrix;
global mg_param;
mg_param.nu1 = 1;            % number of pre smooths
mg_param.nu2 = 1;            % number of post smooths
mg_param.minlevel = 1;       % level of coarsest grid
mg_param.bx  = 5;            % block size x           
mg_param.by  = 5;            % block size y
mg_param.sx  = 4;            % skip size x
mg_param.sy  = 4;            % skip size y
cycleType = 'f';             % MG cycle type.
                             % the other choice is 'v'

% rectangular grids with cell-centered unknowns.
h    = 1/n;
x1d  = [0:h:1];
x1dc = (x1d(1:end-1)+x1d(2:end))/2;
[X,Y] = meshgrid(x1dc, x1dc);

% some RHS function of the test in page 64 of A Multigrid Tutorial
X2 = X.*X;
Y2 = Y.*Y;
RHS = (12*X2-2).*Y2*(Y2-1) + (12*Y2-2).*X2*(X2-1);
rhs = RHS(:);

% the depth of v-cycle.
q=log2(n);
% scale the RHS with h^2 and assign h=1 for the Laplacian operator.
rhs = h*h*rhs;
% assemble the LaplacianOperator
L2d = LaplacianOperator(n, 2, 1);
% the exact solution we expect.
sln = L2d\rhs;
% set the gobal matrices and initialize the MG 
mg_matrix = mg_matrix_init(L2d, q, 'p1');

% use zero as the initial guess.
uS = zeros(n*n,1); 
% the initial residual
res0 = norm(rhs,inf);
% the initial error
err0 = norm(sln,inf);
for jj=1:maxIt
    if cycleType=='v'
%        uS = vcycle_block(uS, rhs, q);
        uS = vcycle(uS, rhs, q);
    elseif cycleType=='f'
%        uS = fcycle_block(uS, rhs, q);
        uS = fcycle(uS, rhs, q);
    else
        error('unknown cycle type!');
    end

    tmp = rhs-L2d*uS;
    res = norm(tmp,inf);
    resRel = res/res0;
    errRel = norm(uS-sln, inf)/err0;
    info = ['  iter #', int2str(jj), ': max-norm of Residual=', ...
        num2str(resRel), '; max-norm of error=', num2str(errRel)];
    disp(info);
%    if (errRel<TOL)
    if (resRel<TOL)
        break
    end
end
nIters = jj;

end
